#################################################################################################################
# Data for the replication of RADIO PUBLIC SERVICE ANNOUNCEMENTS AND VOTER PARTICIPATION AMONG NATIVE AMERICANS: 
#   EVIDENCE FROM TWO FIELD EXPERIMENTS
# Eline A. de Rooij and Donald P. Green
# Political Behavior (POBE)
# Date: July 8, 2016
#################################################################################################################

# install R packages

#install.packages("foreign")
#install.packages("readstata13")

# load R packages

library(foreign)
library(readstata13)

# set working directory   

setwd("") # set working directory


rm(list=ls()) # clear objects in memory   

# read in data frame

subsample <- read.dta13("All stations 2008 and 2010.dta")


##############################
# TABLE 2 MANUSCRIPT
# 2008 and 2010 combined
# ITT, one-tailed p-value 
# and 95%-confidence interval
##############################

outcome <- subsample$outcome
treat <- subsample$treatment
pi <- subsample$pi
npi <- subsample$npi
Treg <- subsample$Treg
total_regna <- subsample$total_regna
state_yr <- subsample$state_yr
exp_round <- subsample$exp_round

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(Treg)
ITT # 0.01860052

treatassn <- function() unlist(sapply(unique(state_yr),function(x) sample(treat[state_yr==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
    ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(Treg)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(Treg)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT)            # 0.14399
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -0.01279620  0.05192591

# widen the confidence interval to reflect the loss of degrees of freedom
inflation_factor <- sqrt( (nrow(as.matrix(Y))-1) / (nrow(as.matrix(Y))-2))
lower <- quantile(ITThyp,0.025)
upper <- quantile(ITThyp,0.975)
center <- (upper + lower)/2 
adj_lower <- center - inflation_factor*(upper-lower)/2
adj_upper <- center + inflation_factor*(upper-lower)/2

c(adj_lower,adj_upper)  # new CI is -0.01299056  0.05212027 

# alternative methods for estimating the standard error of the ITT
# method 1: assess the dispersion of the sampling distribution of the implied schedule of potential outcomes
sd(ITThyp)                      # 0.01691013

# method 2: back out the SE assuming a normal approximation to the sampling distribution
implied_SE <- (adj_upper - adj_lower)/(2*1.96)
implied_SE                      # 0.01660991


##############################
# TABLE 2 MANUSCRIPT
# 2008
# ITT, one-tailed p-value 
# and 95%-confidence interval
##############################

subsample_2008<-subset(subsample, exp_round < 2010)

outcome <- subsample_2008$outcome
treat <- subsample_2008$treatment
pi <- subsample_2008$pi
npi <- subsample_2008$npi
Treg08 <- subsample_2008$Treg08
total_regna <- subsample_2008$total_regna
state_id <- subsample_2008$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(Treg08)
ITT #0.01681381

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(Treg08)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(Treg08)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT)            #  0.2445
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) #  -0.02220010  0.05649373

# widen the confidence interval to reflect the loss of degrees of freedom
inflation_factor <- sqrt( (nrow(as.matrix(Y))-1) / (nrow(as.matrix(Y))-2))
lower <- quantile(ITThyp,0.025)
upper <- quantile(ITThyp,0.975)
center <- (upper + lower)/2 
adj_lower <- center - inflation_factor*(upper-lower)/2
adj_upper <- center + inflation_factor*(upper-lower)/2

c(adj_lower,adj_upper) # new CI is -0.02275824  0.05705187


# alternative methods for estimating the standard error of the ITT
# method 1: assess the dispersion of the sampling distribution of the implied schedule of potential outcomes
sd(ITThyp)                      # 0.02040339
# method 2: back out the SE assuming a normal approximation to the sampling distribution
implied_SE <- (adj_upper - adj_lower)/(2*1.96)
implied_SE                      #  0.02035972


##############################
# TABLE 2 MANUSCRIPT
# 2010
# ITT, one-tailed p-value 
# and 95%-confidence interval
##############################

subsample_2010<-subset(subsample, exp_round > 2008)

outcome <- subsample_2010$outcome
treat <- subsample_2010$treatment
pi <- subsample_2010$pi
npi <- subsample_2010$npi
Treg10 <- subsample_2010$Treg10
total_regna <- subsample_2010$total_regna
state_id <- subsample_2010$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(Treg10)
ITT #0.01986812

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(Treg10)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(Treg10)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) #  0.21841
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -0.02491467  0.06892173

# widen the confidence interval to reflect the loss of degrees of freedom
inflation_factor <- sqrt( (nrow(as.matrix(Y))-1) / (nrow(as.matrix(Y))-2))
lower <- quantile(ITThyp,0.025)
upper <- quantile(ITThyp,0.975)
center <- (upper + lower)/2 
adj_lower <- center - inflation_factor*(upper-lower)/2
adj_upper <- center + inflation_factor*(upper-lower)/2

c(adj_lower,adj_upper)  # new CI is -0.02542190  0.06942897

# alternative methods for estimating the standard error of the ITT
# method 1: assess the dispersion of the sampling distribution of the implied schedule of potential outcomes
sd(ITThyp)                      #  0.02521854
# method 2: back out the SE assuming a normal approximation to the sampling distribution
implied_SE <- (adj_upper - adj_lower)/(2*1.96)
implied_SE                      #  0.02419665



#################################################################################################################
# APPENDIX TABLES
#################################################################################################################

# Note: for the Appendix, we do not widen the confidence interval to reflect the loss of degrees of freedom

##############################
# TABLE A.13 APPENDIX
# 2008 PER STATE
# ITT, one-tailed p-value 
# and 95%-confidence interval
##############################

### AK 2008

AK08<-subset(subsample, state_yr=='AK08')

outcome <- AK08$outcome
treat <- AK08$treatment
pi <- AK08$pi
npi <- AK08$npi
TregAK08 <- AK08$Treg_state08
total_regna <- AK08$total_regna
state_id <- AK08$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregAK08)
ITT # 0.04873096

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(TregAK08)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(TregAK08)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.11326
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # 0.01940787 0.07805406


### AZ 2008

AZ08<-subset(subsample, state_yr=='AZ08')

outcome <- AZ08$outcome
treat <- AZ08$treatment
pi <- AZ08$pi
npi <- AZ08$npi
TregAZ08 <- AZ08$Treg_state08
total_regna <- AZ08$total_regna
state_id <- AZ08$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregAZ08)
ITT # 0.02377179

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(TregAZ08)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(TregAZ08)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.3326
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -0.02025512  0.06779870 


### CA 2008

CA08<-subset(subsample, state_yr=='CA08')

outcome <- CA08$outcome
treat <- CA08$treatment
pi <- CA08$pi
npi <- CA08$npi
TregCA08 <- CA08$Treg_state08

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregCA08)
ITT # 0.06081548


### CO 2008

CO08<-subset(subsample, state_yr=='CO08')

outcome <- CO08$outcome
treat <- CO08$treatment
pi <- CO08$pi
npi <- CO08$npi
TregCO08 <- CO08$Treg_state08

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregCO08)
ITT # -0.3838384


### MN 2008

MN08<-subset(subsample, state_yr=='MN08')

outcome <- MN08$outcome
treat <- MN08$treatment
pi <- MN08$pi
npi <- MN08$npi
TregMN08 <- MN08$Treg_state08
total_regna <- MN08$total_regna
state_id <- MN08$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregMN08)
ITT # -0.008077545

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(TregMN08)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(TregMN08)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.59947
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -0.10370205  0.08735666 


### NM 2008

NM08<-subset(subsample, state_yr=='NM08')

outcome <- NM08$outcome
treat <- NM08$treatment
pi <- NM08$pi
npi <- NM08$npi
TregNM08 <- NM08$Treg_state08
total_regna <- NM08$total_regna
state_id <- NM08$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregNM08)
ITT # 0.05488969

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(TregNM08)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(TregNM08)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.33252
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # 0.02451732 0.08526205


### OK 2008

OK08<-subset(subsample, state_yr=='OK08')

outcome <- OK08$outcome
treat <- OK08$treatment
pi <- OK08$pi
npi <- OK08$npi
TregOK08 <- OK08$Treg_state08

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregOK08)
ITT # -0.1888412


### OR 2008

OR08<-subset(subsample, state_yr=='OR08')

outcome <- OR08$outcome
treat <- OR08$treatment
pi <- OR08$pi
npi <- OR08$npi
TregOR08 <- OR08$Treg_state08

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregOR08)
ITT # 0.04426003


### SD 2008

SD08<-subset(subsample, state_yr=='SD08')

outcome <- SD08$outcome
treat <- SD08$treatment
pi <- SD08$pi
npi <- SD08$npi
TregSD08 <- SD08$Treg_state08

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregSD08)
ITT # -0.1256788


### WA 2008

WA08<-subset(subsample, state_yr=='WA08')

outcome <- WA08$outcome
treat <- WA08$treatment
pi <- WA08$pi
npi <- WA08$npi
TregWA08 <- WA08$Treg_state08
total_regna <- WA08$total_regna
state_id <- WA08$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregWA08)
ITT # -0.06845637

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(TregWA08)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(TregWA08)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.66618
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -0.1767401  0.0398274


### WI 2008

WI08<-subset(subsample, state_yr=='WI08')

outcome <- WI08$outcome
treat <- WI08$treatment
pi <- WI08$pi
npi <- WI08$npi
TregWI08 <- WI08$Treg_state08

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregWI08)
ITT # 0.0989011


### WY 2008

WY08<-subset(subsample, state_yr=='WY08')

outcome <- WY08$outcome
treat <- WY08$treatment
pi <- WY08$pi
npi <- WY08$npi
TregWY08 <- WY08$Treg_state08

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregWY08)
ITT # 0.2828283


### All 2008 - same as TABLE 2 MANUSCRIPT


##############################
# TABLE A.14 APPENDIX
# 2010 PER STATE
# ITT, one-tailed p-value 
# and 95%-confidence interval
##############################

### AK 2010

AK10<-subset(subsample, state_yr=='AK10')

outcome <- AK10$outcome
treat <- AK10$treatment
pi <- AK10$pi
npi <- AK10$npi
TregAK10 <- AK10$Treg_state10
total_regna <- AK10$total_regna
state_id <- AK10$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregAK10)
ITT # 0.02689618

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(TregAK10)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(TregAK10)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.31314
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -0.06012765  0.11392001 


### AZ 2010

AZ10<-subset(subsample, state_yr=='AZ10')

outcome <- AZ10$outcome
treat <- AZ10$treatment
pi <- AZ10$pi
npi <- AZ10$npi
TregAZ10 <- AZ10$Treg_state10
total_regna <- AZ10$total_regna
state_id <- AZ10$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregAZ10)
ITT # 0.04825531

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(TregAZ10)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(TregAZ10)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.16642
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # 0.02859312 0.06791750 


### CO 2010

CO10<-subset(subsample, state_yr=='CO10')

outcome <- CO10$outcome
treat <- CO10$treatment
pi <- CO10$pi
npi <- CO10$npi
TregCO10 <- CO10$Treg_state10
total_regna <- CO10$total_regna
state_id <- CO10$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregCO10)
ITT # -0.03161592

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(TregCO10)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(TregCO10)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.66589
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -0.065465447  0.002233602


### MN 2010

MN10<-subset(subsample, state_yr=='MN10')

outcome <- MN10$outcome
treat <- MN10$treatment
pi <- MN10$pi
npi <- MN10$npi
TregMN10 <- MN10$Treg_state10
total_regna <- MN10$total_regna
state_id <- MN10$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregMN10)
ITT # 0.1189624

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(TregMN10)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(TregMN10)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.25724
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # 0.03508566 0.19215654


### MT 2010

MT10<-subset(subsample, state_yr=='MT10')

outcome <- MT10$outcome
treat <- MT10$treatment
pi <- MT10$pi
npi <- MT10$npi
TregMT10 <- MT10$Treg_state10

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregMT10)
ITT # -0.01846509


### NM 2010

NM10<-subset(subsample, state_yr=='NM10')

outcome <- NM10$outcome
treat <- NM10$treatment
pi <- NM10$pi
npi <- NM10$npi
TregNM10 <- NM10$Treg_state10
total_regna <- NM10$total_regna
state_id <- NM10$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregNM10)
ITT # -0.07217594

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(TregNM10)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(TregNM10)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 1
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -0.1067309 -0.0376210 


### ND 2010

ND10<-subset(subsample, state_yr=='ND10')

outcome <- ND10$outcome
treat <- ND10$treatment
pi <- ND10$pi
npi <- ND10$npi
TregND10 <- ND10$Treg_state10
total_regna <- ND10$total_regna
state_id <- ND10$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregND10)
ITT # 0.2693979

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(TregND10)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(TregND10)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.16642
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -0.002855138  0.541650917


### OK 2010

OK10<-subset(subsample, state_yr=='OK10')

outcome <- OK10$outcome
treat <- OK10$treatment
pi <- OK10$pi
npi <- OK10$npi
TregOK10 <- OK10$Treg_state10

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregOK10)
ITT # 0.1562664


### OR 2010

OR10<-subset(subsample, state_yr=='OR10')

outcome <- OR10$outcome
treat <- OR10$treatment
pi <- OR10$pi
npi <- OR10$npi
TregOR10 <- OR10$Treg_state10

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregOR10)
ITT # -0.1148615


### SD 2010

SD10<-subset(subsample, state_yr=='SD10')

outcome <- SD10$outcome
treat <- SD10$treatment
pi <- SD10$pi
npi <- SD10$npi
TregSD10 <- SD10$Treg_state10
total_regna <- SD10$total_regna
state_id <- SD10$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregSD10)
ITT # 0.03246393

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(TregSD10)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(TregSD10)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.33246
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # 0.03168562 0.03324223


### WA 2010

WA10<-subset(subsample, state_yr=='WA10')

outcome <- WA10$outcome
treat <- WA10$treatment
pi <- WA10$pi
npi <- WA10$npi
TregWA10 <- WA10$Treg_state10
total_regna <- WA10$total_regna
state_id <- WA10$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregWA10)
ITT # -0.02589112

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(TregWA10)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(TregWA10)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.64862
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -0.13499109  0.08320886


### WI 2010

WI10<-subset(subsample, state_yr=='WI10')

outcome <- WI10$outcome
treat <- WI10$treatment
pi <- WI10$pi
npi <- WI10$npi
TregWI10 <- WI10$Treg_state10

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregWI10)
ITT # -0.0471597


### WY 2010

WY10<-subset(subsample, state_yr=='WY10')

outcome <- WY10$outcome
treat <- WY10$treatment
pi <- WY10$pi
npi <- WY10$npi
TregWY10 <- WY10$Treg_state10

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(TregWY10)
ITT # 0.3459038


### All 2010 - same as TABLE 2 MANUSCRIPT


##############################
# TABLE A.15 APPENDIX
# NON-NATIVE STATIONS 
# ITT, one-tailed p-value 
# and 95%-confidence interval
##############################

rm(list=ls()) # clear objects in memory   

# read in data frame

Nonnat <- read.dta13("Nonnative stations 2008 and 2010.dta")


# 2008 and 2010 combined

outcome <- Nonnat$outcome
treat <- Nonnat$treatment
pi <- Nonnat$pi
npi <- Nonnat$npi
Treg_nnstation <- Nonnat$Treg_nnstation
total_regna <- Nonnat$total_regna
state_yr <- Nonnat$state_yr
exp_round <- Nonnat$exp_round

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(Treg_nnstation)
ITT # 0.02842157

treatassn <- function() unlist(sapply(unique(state_yr),function(x) sample(treat[state_yr==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y<- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(Treg_nnstation)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(Treg_nnstation)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull) 
mean(ITTnull >= ITT) # 0.23515
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -0.01873477  0.07853936 


# 2008 

Nonnat_2008<-subset(Nonnat, exp_round < 2010)

outcome <- Nonnat_2008$outcome
treat <- Nonnat_2008$treatment
pi <- Nonnat_2008$pi
npi <- Nonnat_2008$npi
Treg_nnstation08 <- Nonnat_2008$Treg_nnstation08
total_regna <- Nonnat_2008$total_regna
state_id <- Nonnat_2008$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(Treg_nnstation08)
ITT # 0.01789113

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y<- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(Treg_nnstation08)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(Treg_nnstation08)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull) 
mean(ITTnull >= ITT) # 0.40771
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -0.03831911  0.08598124 


# 2010

Nonnat_2010<-subset(Nonnat, exp_round > 2008)

outcome <- Nonnat_2010$outcome
treat <- Nonnat_2010$treatment
pi <- Nonnat_2010$pi
npi <- Nonnat_2010$npi
Treg_nnstation10 <- Nonnat_2010$Treg_nnstation10
total_regna <- Nonnat_2010$total_regna
state_id <- Nonnat_2010$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(Treg_nnstation10)
ITT # 0.0523876

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y<- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(Treg_nnstation10)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(Treg_nnstation10)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull) 
mean(ITTnull >= ITT) # 0.12362
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -0.0294617  0.1341139


##############################
# TABLE A.15 APPENDIX
# NATIVE STATIONS 
# ITT, one-tailed p-value 
# and 95%-confidence interval
##############################

rm(list=ls()) # clear objects in memory   

# read in data frame

Nat <- read.dta13("Native stations 2008 and 2010.dta")

# 2008 and 2010 combined 

outcome <- Nat$outcome
treat <- Nat$treatment
pi <- Nat$pi
npi <- Nat$npi
Treg_nastation <- Nat$Treg_nastation
total_regna <- Nat$total_regna
state_yr <- Nat$state_yr
exp_round <- Nat$exp_round

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(Treg_nastation)
ITT # 0.03432274

treatassn <- function() unlist(sapply(unique(state_yr),function(x) sample(treat[state_yr==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y<- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(Treg_nastation)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(Treg_nastation)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull) 
mean(ITTnull >= ITT) # 0.02394
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # 0.008519487 0.062146765


# 2008 

Nat_2008<-subset(Nat, exp_round < 2010)

outcome <- Nat_2008$outcome
treat <- Nat_2008$treatment
pi <- Nat_2008$pi
npi <- Nat_2008$npi
Treg_nastation08 <- Nat_2008$Treg_nastation08
total_regna <- Nat_2008$total_regna
state_id <- Nat_2008$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(Treg_nastation08)
ITT # 0.02488196

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y<- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(Treg_nastation08)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(Treg_nastation08)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull) 
mean(ITTnull >= ITT) #  0.24903
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -0.01843802  0.06820194


# 2010

Nat_2010<-subset(Nat, exp_round > 2008)

outcome <- Nat_2010$outcome
treat <- Nat_2010$treatment
pi <- Nat_2010$pi
npi <- Nat_2010$npi
Treg_nastation10 <- Nat_2010$Treg_nastation10
total_regna <- Nat_2010$total_regna
state_id <- Nat_2010$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/mean(Treg_nastation10)
ITT # 0.03803227

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y<- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT*total_regna[treat==1]
Y1h[treat == 0] <- Y[treat==0] + ITT*total_regna[treat==0]

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/mean(Treg_nastation10)
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/mean(Treg_nastation10)
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull) 
mean(ITTnull >= ITT) # 0.062
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # 0.007746789 0.066563208


##############################
# TABLE A.16 APPENDIX
# ALL STATIONS 
# ITT, one-tailed p-value 
# and 95%-confidence interval
##############################

rm(list=ls()) # clear objects in memory   

# read in data frame

subsample <- read.dta13("All stations 2008 and 2010.dta")


# 2008 and 2010 combined

outcome <- subsample$outcome
treat <- subsample$treatment
pi <- subsample$pi
npi <- subsample$npi
Treg <- subsample$Treg
total_regna <- subsample$total_regna
state_yr <- subsample$state_yr
exp_round <- subsample$exp_round

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/85
ITT # 23.77059

treatassn <- function() unlist(sapply(unique(state_yr),function(x) sample(treat[state_yr==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT
Y1h[treat == 0] <- Y[treat==0] + ITT

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/85
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/85
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.14399
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -16.66429  66.59586 


# 2008

subsample_2008<-subset(subsample, exp_round < 2010)

outcome <- subsample_2008$outcome
treat <- subsample_2008$treatment
pi <- subsample_2008$pi
npi <- subsample_2008$npi
Treg08 <- subsample_2008$Treg08
total_regna <- subsample_2008$total_regna
state_id <- subsample_2008$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/37
ITT # 20.48649

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT
Y1h[treat == 0] <- Y[treat==0] + ITT

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/37
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/37
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.2445
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -28.89204  74.18030


## 2010

subsample_2010<-subset(subsample, exp_round > 2008)

outcome <- subsample_2010$outcome
treat <- subsample_2010$treatment
pi <- subsample_2010$pi
npi <- subsample_2010$npi
Treg10 <- subsample_2010$Treg10
total_regna <- subsample_2010$total_regna
state_id <- subsample_2010$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/48
ITT # 26.30208

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y <- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT
Y1h[treat == 0] <- Y[treat==0] + ITT

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/48
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/48
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull)
mean(ITTnull >= ITT) # 0.21841
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -31.63660  89.58167


##############################
# TABLE A.16 APPENDIX
# NONNATIVE STATIONS 
# ITT, one-tailed p-value 
# and 95%-confidence interval
##############################

rm(list=ls()) # clear objects in memory   

# read in data frame

Nonnat <- read.dta13("Nonnative stations 2008 and 2010.dta")


# 2008 and 2010 combined

outcome <- Nonnat$outcome
treat <- Nonnat$treatment
pi <- Nonnat$pi
npi <- Nonnat$npi
Treg_nnstation <- Nonnat$Treg_nnstation
total_regna <- Nonnat$total_regna
state_yr <- Nonnat$state_yr
exp_round <- Nonnat$exp_round

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/41
ITT #  21.14634

treatassn <- function() unlist(sapply(unique(state_yr),function(x) sample(treat[state_yr==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y<- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT
Y1h[treat == 0] <- Y[treat==0] + ITT

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/41
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/41
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull) 
mean(ITTnull >= ITT) # 0.23515
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -22.97829  66.99882


# 2008 

Nonnat_2008<-subset(Nonnat, exp_round < 2010)

outcome <- Nonnat_2008$outcome
treat <- Nonnat_2008$treatment
pi <- Nonnat_2008$pi
npi <- Nonnat_2008$npi
Treg_nnstation08 <- Nonnat_2008$Treg_nnstation08
total_regna <- Nonnat_2008$total_regna
state_id <- Nonnat_2008$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/19
ITT # 19.95614

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y<- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT
Y1h[treat == 0] <- Y[treat==0] + ITT

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/19
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/19
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull) 
mean(ITTnull >= ITT) # 0.40771
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -57.42805 108.63612


# 2010

Nonnat_2010<-subset(Nonnat, exp_round > 2008)

outcome <- Nonnat_2010$outcome
treat <- Nonnat_2010$treatment
pi <- Nonnat_2010$pi
npi <- Nonnat_2010$npi
Treg_nnstation10 <- Nonnat_2010$Treg_nnstation10
total_regna <- Nonnat_2010$total_regna
state_id <- Nonnat_2010$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/22
ITT # 22.17424

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y<- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT
Y1h[treat == 0] <- Y[treat==0] + ITT

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/22
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/22
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull) 
mean(ITTnull >= ITT) #  0.12362
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -12.70018  57.04362


##############################
# TABLE A.16 APPENDIX
# NATIVE STATIONS 
# ITT, one-tailed p-value 
# and 95%-confidence interval
##############################

rm(list=ls()) # clear objects in memory   

# read in data frame

Nat <- read.dta13("Native stations 2008 and 2010.dta")

# 2008 and 2010 combined 

outcome <- Nat$outcome
treat <- Nat$treatment
pi <- Nat$pi
npi <- Nat$npi
Treg_nastation <- Nat$Treg_nastation
total_regna <- Nat$total_regna
state_yr <- Nat$state_yr
exp_round <- Nat$exp_round

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/19
ITT # 85.44737

treatassn <- function() unlist(sapply(unique(state_yr),function(x) sample(treat[state_yr==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y<- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT
Y1h[treat == 0] <- Y[treat==0] + ITT

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/19
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/19
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull) 
mean(ITTnull >= ITT) #  0.02394
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # 8.466758 162.298475


# 2008 

Nat_2008<-subset(Nat, exp_round < 2010)

outcome <- Nat_2008$outcome
treat <- Nat_2008$treatment
pi <- Nat_2008$pi
npi <- Nat_2008$npi
Treg_nastation08 <- Nat_2008$Treg_nastation08
total_regna <- Nat_2008$total_regna
state_id <- Nat_2008$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/6
ITT # 55.33333

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y<- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT
Y1h[treat == 0] <- Y[treat==0] + ITT

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/6
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/6
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull) 
mean(ITTnull >= ITT) #  0.24903
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -20.44444 131.11111


# 2010

Nat_2010<-subset(Nat, exp_round > 2008)

outcome <- Nat_2010$outcome
treat <- Nat_2010$treatment
pi <- Nat_2010$pi
npi <- Nat_2010$npi
Treg_nastation10 <- Nat_2010$Treg_nastation10
total_regna <- Nat_2010$total_regna
state_id <- Nat_2010$state_id

ITT <- (sum(outcome[treat==1]/pi[treat==1]) - sum(outcome[treat==0]/npi[treat==0]))/13
ITT #  99.34615

treatassn <- function() unlist(sapply(unique(state_id),function(x) sample(treat[state_id==x])))

set.seed(1234567)

numiter <- 100000

ITTnull <- ITThyp <- rep(NA,numiter)

Y0h <- Y1h <- Y<- outcome
Y0h[treat == 1] <- Y[treat==1] - ITT
Y1h[treat == 0] <- Y[treat==0] + ITT

for (iter in 1:numiter) {
  treatri <- treatassn()
  ITTnull[iter] <- (sum(Y[treatri==1]/pi[treatri==1]) - sum(Y[treatri==0]/npi[treatri==0]))/13
  ITThyp[iter] <- (sum(Y1h[treatri==1]/pi[treatri==1]) - sum(Y0h[treatri==0]/npi[treatri==0]))/13
  if(iter %% 10000 == 1) cat(iter,"")
}

summary(ITTnull) 
mean(ITTnull >= ITT) #  0.062
summary(ITThyp)
quantile(ITThyp,c(0.025,0.975)) # -13.86834 211.83284 

